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Abstract 

This report presents a computational technique for the rapid, efficient calculation of 
fields from transient acoustic sources in linear, isotropic media. The source velocity is 
separable in space and time. The method uses a spatial impulse response method based on 
linear systems concepts to express the output in terms of the Green’s function of propagation 
equation and the boundary conditions. The output is expressed as the inverse spatial 
transform of the product of the transform of the spatial excitation and a time-varying 
spatial filter that represents propagation. The calculation technique has been implemented 
in MATLAB and sample cases are presented for the circular and square piston, as well as a 
Gaussian- and Bessel-weighted spatial excitation. Code for the MATLAB implementation 
is provided. 
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Chapter 1 

Introduction 



1.1 Introduction 

The propagation characteristics of continuously radiated monochromatic ultrasonic sources 
are well solved through application of the angular spectrum technique [l] or Fresnel in- 
tegrals. More frequently, however, acoustic imaging, tissue characterization, and physical 
acoustics applications tend to use pulsed sound. The propagation of pulsed ultrasound with 
arbitrary temporal and spatial components is not understood to the same degree. We want 
to develop a reliable, easy method of diffraction prediction. This report describes an ap- 
proach based on linear systems theory and the Fourier transform. The goal was to achieve a 
readily usable method of predicting pulsed wave diffraction in a time-efficient and accurate 
manner in order to examine the wave diffraction. Earlier efforts had used FORTRAN code 
to implement the propagation model on both a mainframe computer [2]- [6] and a personal 
computer [7]. Merrill had attempted to use MATLAB (a commercial matrix manipulation 
software package) to perform the calculations on a PC but was thwarted by memory re- 
strictions of the earlier versions of this software. More recent versions allow larger matrices 
and arrays to be manipulated (subject only to the computer memory available) Thesis work 
performed by Upton [8] and Reid [9] applied this updated software to the problems of op- 
tical propagation and acoustic propagation, respectively. Much of the foundation for the 
work reported here was done by these officer-students. 

The basic method of the spatial impulse response was introduced by Stepanishen [10]- 
[13], reviewed by Harris [14], and modified by Guyomar and Powers [5]. The approach of 
Guyomar and Powers differed by the use of linear systems theory. Linear systems theory 
revealed the importance of the total impulse response and its equivalence to the Green’s 
function. Furthermore, the spatial impulse response functions are found in the spatial 
transform (Fourier) domain. Working in the transform domain allowed propagation of the 
wave to be viewed as a time-varying spatial filter applied to the spatial spectrum of the 
input excitation. The advantage of this method is that it provides the diffracted field from 
an insonifying wave with arbitrary temporal and spatial dependence in a computationally 
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efficient form. By use of the Fourier transform, an efficient computer implementation of 
this technique using FFT routines is possible. 

The desired benefit of a fast, time-efficient computer implementation to calculate the 
acoustic potential or pressure is to aid in ultrasonic transducer design for medical, acoustic 
imaging, and mine warfare applications. With a knowledge of wave diffraction phenomenon 
a diffracted wave reflected from an unknown object can be used to provide information 
about the object. This type of system must be portable as well as time-efficient, which is 
very achievable given the trends in computer technology. Computers have become faster and 
have increased memory capacity while their size has decreased. Other benefits are derived 
from the use of the matrix manipulation program, MATLAB, and the ability to expand 
this implementation to cases involving lossy media. Because MATLAB is readily available 
on the commercial market, it requires no special equipment for computer implementation. 



1.2 Report Overview 

Chapter II consists of the problem description, including the source-to-observation plane 
geometry, a discussion on the linear systems approach, the mathematical development, and 
an overview of MATLAB and AXUM. Chapter III consists of a discussion of the two pro- 
gram modules, the acoustic filter module, ACJFIL.M, and the acoustic propagation module, 
AC-PROP.M. Chapter IV starts with the set of defining parameters. These parameters are 
then used for an explanation of the program’s verification and an investigation of other 
input excitation functions. Chapter V contains some examples of the numerical output for 
various source excitation conditions. 

Following a summary in Chapter VI, Appendices A and C give detailed explanations 
of the MATLAB routines that were used to produce the numerical results, AC-FIL.M and 
ACJPROP.M. Appendix B gives examples of the propagation filters generated by the code 
in Appendix A. The source code of the various geometric excitation functions is given in 
Appendix D. Representative output presentations using a scientific visualization program, 
called Spyglass Dicer, are found in Appendix E. 
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Chapter 2 

Problem Description 



Before assembling a computer implementation, we must first understand the problem. An 
explanation of the problem follows, beginning with the description of the geometry in the 
first section. Section two continues the explanation into the linear systems approach. The 
third section proceeds through the mathematical development of the problem and ties in 
the Fourier transform. The theory presented in the first three sections was derived from 
the works of Guyomar and Powers [3, 4, 5, 6]. The final section gives an overview of the 
software tools used for generating the excitation functions, the propagation fields, and the 
output graphics. 

2.1 Geometry 

The problem geometry is shown in Fig. 2.1. The acoustic velocity potential 0(z,y,z) is 
to be calculated at an arbitrary point in the positive- z half-space given the z-directed 
velocity in a portion of the source plane. The source’s z-directed velocity is spatially and 
temporally arbitrary and is rigidly baffled (i.e., equal to zero) in the region outside the 
source. Furthermore, it is assumed that the spatial and temporal components of the z- 
velocity are separable, having the form at the input plane 

v*(*o,yo,(M) = T(t)s(z 0 ,i/o). (2.1) 

A linear, homogeneous, and lossless medium is assumed to be present between the source 
and observation point. (The extension to lossy propagation can be found in Ref. [15].) 

2.2 Linear Systems Approach 

Linear systems solutions are applied to systems that are linear and time-invariant. A 
linear systems solution approach to this problem is possible because propagation in a linear 
homogeneous media is a linear, space-invariant process [3]. In linear systems the impulse 
response is the response of the system to an impulsive input. The total impulse response 
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/i(x,y,z,f) of a system is produced by an input of the form <5(x,j/,*) = 6{x,y)6{t)] this is 
shown in Fig. 2.2a. Figure 2.2b shows the spatial impulse response p(x,j/,z,t), defined as 
the response to an excitation of the form s(x,y)b(t). Note that an arbitrary spatial input 
has been substituted for the impulsive spatial input. Recall from linear systems theory that 
the solution for an arbitrary input (spatial or temporal) is the convolution of the input with 
the system’s total impulse response; therefore, the spatial impulse response in this case is 
given by 



p(x,y,z,t ) = s{x,y)‘ x m y h(x,y,z,t), (2.2) 

where * indicates convolution with respect to the variable shown. 

The double convolution in Eq. 2.2 is not easily computer implemented. The Fourier 
transform furnishes a convenient method to resolve this dilemma by using the property 
that convolution in the spatial domain becomes multiplication in the transform domain, 
i.e., 

p(fxify> z tt) = ^(/r> fy )M/r> fy > z , t) , (2-3) 

where the tilde notation indicates the Fourier transform of the function (in this case, the 
two-dimensional spatial Fourier transform). This relation is shown in Fig. 2.2c. 

The general solution is shown in Fig. 2. 2d, where 

4>(x,y,z,t) = s{x,y)T(t)l m y m t h(x,y,z,t). (2.4) 
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,y 4h & boundary ' > 

conditions 
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Figure 2.2: Block diagram of (a) the total impulse response h(x,y,t)), (b) the spatial 
impulse response p(x,y,z,t) in the space-time domain, (c) the spatial impulse response 
P(fx,fy,t ) in the spatial frequency-time domain, and (d) the general solution p(x,y,z,t). 
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Substituting Eq. 2.2 into Eq. 2.4 gives 

<t>{x,y,z,t) = T(t);p(x,y,z,t). (2.5) 

It follows from Fig. 2.2c and Eqs. 2.2 and 2.5 that the key to finding the general solution 
is finding the spatial impulse response p(x,y,z,t) which is, in turn, dependent on the total 
impulse response of the system /i(x, y, z, *). This total impulse response of the system is 
the propagation field that results from an impulsive source, as in Fig 2.2a, that solves the 
wave equation and satisfies the boundary conditions. The solution to the wave equation 
satisfying the boundary conditions is commonly known as the Green’s function ; hence, the 
total impulse response is simply the Green’s function. Therefore, once the Green’s function 
is known, the total impulse response is also known, and the solution becomes a triple 
convolution between an excitation source which is spatially and temporally arbitrary (and 
assumed to be known), and the systems’ total impulse response or Green’s function. The 
two spatial convolutions are easily implemented in the spatial transform domain; the time 
convolution can be implemented in the temporal transform domain, if desired. We have 
found it easier to perform the time convolution directly in the time domain. 



2.3 Mathematical Development for Acoustic Propagation 



The wave equation for lossless media is 

, 1 d 2 <t> 

- °- (2 - 6 > 

The general solution of Fig. 2.2c gives the result in terms of the acoustic potential <p which 
must be found from a z - velocity input. We relate the acoustic velocity and acoustic potential 
with the following relationship; 



v(x,y,z,t) 



-V(t>(x,y, z,t) . 



(2.7) 



Hence, the 2 -velocity component is given by 

v z {x,y,z,t) = 



dz 



( 2 . 8 ) 



Since the wave equation is in terms of c (the acoustic velocity in the media) and time t, the 
partial derivative with respect to z must be related to these two parameters. This is done 
by using the fact that a wave traveling in the positive-z direction has an argument of the 
form ct — z, resulting in the relationship 



d_ _ d_ 

dz °dt ' 



(2.9) 



Applying Eq. 2.9 to Eq. 2.8 gives the z-velocity at the input plane (z = 0) as 



v z (x,y,z,t) 



d<t>(x,y,o,t) 
c dt 



(2.10) 
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Equation 2.10 requires the acoustic potential at the input plane. It was assumed in 
Eq. 2.1 that the z-velocity is separable which also implies a separable acoustic potential. 
Such an acoustic potential has the general form 

4>{x,y,0,t) = s(x,p)r(f) (2.11) 

at the z — 0 plane. If Eq. 2.11 is substituted into Eq. 2.10 and the partial derivative is 
carried out, then Eq. 2.10 becomes 

Vz(x,y,z,t) = cs(x,y)T'(t), (2.12) 

where the prime indicates a derivative with respect to the time variable t. A comparison 
of the z - velocities given in Eqs. 2.1 and 2.12 indicates that T(t) is equivalent to cr'(t). 
Equation 2.12 is now the input in Fig. 2.2c. 

As stated earlier, the general solution of Fig. 2.2c is the Green’s function. For the 
standard wave equation (for lossless media) and rigidly baffled boundary conditions, the 
applicable Green’s function is [4] 

h{x,y,z,t) = (2-13) 

where R = \Jx 2 + y 2 + z 1 . Substituting this into Eq. 2.2 provides the spatial impulse 
response 



p(x,y,z,t) = s(x y y) m x ’ h(x y y,z y t) 

= s{x,y) x y j . (2.14) 

Thus, the spatial impulse response in the spatial transfer domain is the product of 
the angular spectrum of the source s and the propagation transfer function p, written 
symbolically as 



P{fx,fy,Z,t) = s(fx,fy,t)h(f x J y ,Z,t). 



(2.15) 



Taking the two-dimensional spatial Fourier transform of the Green’s function in Eq. 2.13 
gives the propagation transfer function 



Syi Z 1 1 ) 



where the relationship 



T 





2,/q H(ct-z), 



Jo (pv/c 2 t 2 - 2 2 ) H{ct - z) 

(et)"- 1 



(2.16) 



(2.17) 
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was applied (with r = y/J + /y). The term H{ct — z) is the Heaviside step function and 
makes the wave causal. 

Equation 2.16 exhibits two important points. The first is that the propagation transfer 
function is a Bessel of the first kind of zero order. Secondly, the propagation transfer 
function can be identified as a time-varying spatial filter having a Bessel shape. The filter 
begins as an all-pass filter early in its time evolution and then becomes more of a low-pass 
filter as time progresses. This increasingly narrow low-pass filter reduces the resolution 
in the receiving plane of the system as time progresses. In our computer programs, the 
propagation spatial filter at various instants of time are produced by the MATLAB file, 
AC-FIL.M. 

Equation 2.16 is valid for an input that is temporally impulsive and spatially arbitrary. 
Taking the inverse two-dimensional spatial transform of Eq. 2.16 would, in this case, give 
a final result since convolution with an impulse is the same function at the time that the 
impulse occurred. To account for a non-impulsive time input component, the convolution 
of Eq. 2.5 must be carried out, resulting in 

<t>(x,y,z,t ) = , (2.18) 

where p is the product of the angular spectrum of the source s and the propagation trans- 
fer function h. Since only temporally impulsive inputs are simulated in the examples that 
follow, the convolution in Eq. 2.18 was not computed for our cases. Our final solution, 
therefore, reduces to taking the inverse 2-D spatial Fourier transform of the spatial im- 
pulse response. The program module, AC_PROP.M, simulates the spatial impulse response 
solution for the input excitation function distribution chosen by the user. 

The program that implements these equations is discussed in Chapter III and detailed 
in Appendix C. Illustrative examples of the time-varying filters and outputs follow in the 
Chapter IV discussion with more examples supplied in Appendices B and E. The tool 
of implementation for the program was MATLAB with graphical assistance provided by 
AXUM. Output data was also represented using a scientific visualization program called 
Spyglass Dicer. An overview of these programs follows in the next section. 



2.4 Software Tools 
2.4.1 MATLAB Overview 

MATLAB is a high-performance, interactive numeric computation software package for 
science and engineering applications produced by The MATHWORKS, Inc. The name 
comes from MATrix LABoratory; hence, the basic data element is a matrix (or array), 
which does not require dimensioning. MATLAB. A major advantage of MATLAB is that 
it uses a “vectorized’ 1 approach to computations, simplifying the programming. Another 
distinct advantage is the availability of preprogrammed functions, such as calculation of the 
two-dimensional FFTs and the Bessel function [16]. 
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In MATLAB there are two type of macro-like files. A script file is used to automate long 
sequences of commands including functions. Arguments are not passed into script files. A 
function file , however, may have arguments passed into it and out of it. Examples of script 
files in this thesis include ACJFIL.M and AC-PROP.M. Examples of functions include the 
input functions, the three dimensional graphing function mesh, and the two “fit” functions 
that realize the Fourier transform [16]. 

The two “fft” functions employed for the Fourier transform are fft2 and fftshift. The 
fft2 function is a two-dimensional fast Fourier transform and is repeatedly used throughout 
our program. MATLAB allows the use of only positive indices in an array. Geometrically, 

this is equivalent to working only in the first quadrant of an £ y plane. Our sources are 

symmetric around the origin and we want to display them using all four quadrants. The 
key to the transformation is the fftshift function. Figure 2.3a. shows the first-quadrant 
representation of a rectangular source centered on the origin. (Note the MATLAB assumes 
implicitly that the function is periodic in the x and y directions. Periodic translation of the 
first-quadrant in the x and y directions will fill in the rest of the source at the origin.) The 
fftshift function rearranges the data so that it is centered in the first quadrant. (Actually 
the data is not quite centered when one uses an even number of sample points, since there are 
not any sample points aligned with the center. The shifted array is offset from the quadrant 
center by one-half of a sample interval in both the x and y directions.) Figure 2.3b represents 
the shifted source. For viewing purposes, the shifted version is easier to understand; for 
computational purposes the unshifted function must be used to avoid phase errors when 
calculating the transform or its inverse. We will use the terms “shifted” for functions 
represented in the centered geometry and “unshifted” for the functions represented in the 
corner geometry. The interested reader is referred the MATLAB User’s Guide [16] for more 
details. 

The MATLAB mesh command can be used to provide three-dimensional plots of the 
computed fields. Alternatively, the input and output data can be stored in ASCII format, 
imported into a more powerful graphics package, called AXUM, and plotted. 

2.4.2 AXUM Overview 

AXUM is an interactive software package for technical graphics and data analysis, produced 
by TriMetrix Inc. AXUM allows easy manipulation of raw data imported from MATLAB 
in ASCII format. (Data may also be imported from several other formats.) Once imported, 
data manipulation can be carried out by AXUM through use of the Transform and Convert 
menu options. Then the data is arranged using the MAT2GRID routine found in the History 
Editor menu, so that the desired surface plot can be generated. The MAT2GR1D routine 
produces three columns of data from the original data array by placing each element of 
acoustic potential data in the z data column with the indices from the array in the x and 
y columns. Once processed, the data is ready for graphing [17]. 

The Graph menu gives various options for controlling the graph’s attributes and general 
aesthetics. Once the desired axes intervals, labels, and titles have been set, the graph can 
be displayed in the Edit Screen window. The Edit Screen window allows changes to be 
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(b) 



Rgme 2.3: Illustration of tie 



center versus corner geometry. 
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made to the graph while it is displayed. (This can be is time-consuming because the graph 
is redrawn after each change.) A useful feature of the Edit Screen window, however, is the 
capability to rotate and tilt the graph interactively so that the preferred aspect is achieved. 
After establishing the desired attributes, the graph can be saved as a graph or as an image. 
The concept of saving the graphs as images is another very useful feature because it allows 
a single graph to be assembled and saved as a graph template on which the other images 
may be overlaid. 

2.4.3 Spyglass Software Overview 

The Spyglass software consists of a set of commercial programs produced by Spyglass, Inc. 
for the purpose of aiding in visualization of scientific data on Macintosh computers. (One 
program, the Transform program, is also commercially available for Unix workstations.) 

The Spyglass Data Utility has the capability to read data from a wide variety of formats 
into the HDF (Hierarchical Data Format) format 1 used by the Spyglass software. In 
particular, for our work, the utility is able to read ASCII data. In our case, we took the 6^ 
ASCII data files produced by propagation model and transferred them to a Macintosh II 
computer. The Spyglass Data Utility was used to create an HDF file, representing a data 
volume of 128x128x64 samples. (The typical data compression was from 16 MB of ASCII 
data to a 6-MB HDF file.) 

Spyglass Dicer is a program that allows the user to interactively look at the 3-D data 
volume in various ways. Planes of any arbitrary orientation can be inserted or rectan- 
gular regions of the data can be observed. (Figures 4.3 and 4.4 on pages 24 and 25 are 
representative of the data views that can be obtained with the Dicer program.) 

The Spyglass Transform program allows the user to pick a plane of the HDF data volume 
(currently, the planes must be oriented perpendicular to one of the axes) and to view it in 
a variety of formats, such as a color intensity display, a gray-scale intensity display, or a 
surface plot. 



*The HDF format is 
platforms. 



becoming a popular means of transferring large data files between computing 
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Chapter 3 

MATLAB Modeling of Equations 



This chapter discusses the two MATLAB files that implement the model, AC-F1L.M and 
AC-PROP.M. Together, these script files form the two modules of the program that simulate 
acoustic wave diffraction via implementation of the concepts and equations presented in the 
previous chapter. The program was written in tw r o modules since the propagation filter 
characteristics of the medium depend only on the locations of the source and receiving 
plane locations and are independent of the source geometry. The propagation filters could 
be computed once and then be reused for a variety of excitation conditions. This modularity 
of approach was fortuitous, since the calculation of the Bessel functions in the filters was 
computationally intensive. (This modular approach to the calculation of the fields is one 
of the major advantages of our method.) A working narrative of AC-FIL.M opens the 
discussion, followed by a working narrative of AC-PROP.M. The excitation functions, or 
input functions, are included under the AC-PROP.M heading. A brief summary of the 
program steps follows in the final section. 



3.1 Acoustic Filter Module 

The script file, AC-FIL.M (ACoustic FILter) computes the time-varying Bessel filters of 
Eq. 2.16, repeated here for convenience as 

M/., /„*,<) = 

= 2 J 0 (p\/c 2 * 2 - z 2 ) H(ct ~ z ) . (3.1) 

Before discussing the coding of Eq. 3.1, the array geometry must be explained. 

The array geometry is set by the size of the array and the sampling frequency. The 
variable base denotes the number of points on a side in the base array creating an base x base 
array. (It is advantageous in calculating an FFT to make base a power of 2; this is not 
necessary, however.) Initially base was given a value of 64 as in previous work [3] — [7]. 
Once the program result was verified, base was increased by a factor of two to 128. Making 



12 



base 



bose 



+1 



Suborroy II 


Suborroy 1 


Subarray HI 


'^''(NO.NO) 
Subarray IV 



bose 



+1 



base 



Array index 

Figure 3.1: Offset geometry of base array matrix. 



base an even number, however, causes the center of symmetry of the array to fall between 
array elements. The center of symmetry of the source A 7 0 coincided with the array element 
at 



A 7 0 



base 

~ 2 ~ 



+ 1 . 



(3.2) 



This results in an offset geometry as shown in Figure 3.1. 

The offset center at (NO, NO) divides the base array into four subarrays, each of different 
sizes. The fact that the subarrays have different sizes is important in the use of symmetry 
because only one subarray is actually computed; the other subarrays are determined by 
symmetry. We calculate the data for the A 7 0 x (A 7 0 — 1) points in Subarray I, shown in 
Fig. 3.1. In addition, data is calculated for an additional column (A 7 0 x 1) located at the 
right side of Subarray I. (The extra column was required to compute all of the values in the 
other three quadrants using symmetry.) Subarrays II, III, and IV of Fig. 3.1 are determined 
from Quadrant I using symmetry with MATLAB’s flip commands. The use of symmetry in 
this manner was employed in both program modules. It should be noted that this assumed 
i- and y-axis symmetry imposes an other restriction on the assumed source geometry; this 
restriction can be removed by performing the computations over the full dimensions of the 
array rather than relying on symmetry to simplify the calculations. 

The number of time samples (or number of time slices), slices, was set at 64. Initially, 
this was to emulate previous work [3] — (7]; however, the time duration and resolution 
proved to be adequate for the follow-on simulations. Although there are 64 time slices, 
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only 61 filters are generated by the MATLAB implementation of Eq. 3.1. The Heaviside 
step function in Eq. 3.1 was simulated by arbitrarily setting a variable Step to three which 
produces three all-zero rows for the first three time slices. The result is that slices - Step, 
or 61, time slices actually require computing. Since the first three time slices are zero, the 
computations start with the fourth time slice at time z/c and proceed to the time value, 
time-max , provided by the user. 

Referring back to Eq. 3.1, it is seen that the argument to the Bessel function, p\/c 2 t 2 -/z 2 , 
is composed of four variables. Of the four, only the time t varies within the program. In 
the MATLAB code, time t is represented by the variable time. As previously stated, filter 
generation does not start until time z/c and is linearly incremented to the maximum time 
of propagation time-max. The source-to-receiver distance z was assigned the value z — 10 
cm in the MATLAB code and c, the acoustic velocity in the medium, was assigned the value 
c = 1500 m/s (velocity in water). The value of z was originally chosen to parallel the pre- 
vious work [3] — [7] and was found to be convenient for Subsequent simulations. The. value 
of time-max was set to 150 /is for the verification phase and to 375 ps for the remainder 
of the simulations. Consequently, time ranged from 66.667 ps to 150 ps (or 375 /is) in 61 
increments. 

The finzil variable in the Bessel argument to be examined is p. In the MATLAB code, p 
was given the name rho and has a maximum value, rho-max , of 200. The value of rho.max 
= 200 was arbitrarily chosen subject to the constraint that it needs to be large enough that 
the field does not extend beyond the edges of the array at its largest width (or else the field 
will be aliased back into the array from the edges). The value of 200 was chosen following 
Merrill [7]. Although rho is not time-varying, it does vary with spatial frequency 

P = v//x 2 + / y 2 - (3-3) 

A vector having NO - 1 points extending from 0 to rho-max was then formed. Then, 
the vector was used in the MATLAB routine meshdom [16] to form two identical matrices 
rhox and rhoy. The meshdom routine takes a given vector and forms two matrices with the 
same spacing in the x direction and y direction. The two matrices, rhox and rhoy represent 
f x and f y in Eq. 3.3 and Fig. 3.2, which shows graphically the construction of rho. In 
Fig. 3.2, the inner scales are in terms of the column and row numbers, respectively. The 
row index runs from 1 to NO ; the column index runs from NO to base. The x and y axes 
can be rescaled to units of spatial frequency ( f x and f y ) with units of cycles/m. This is 
represented in Fig. 3.2 by the outer labeling. 

The combination of rho and the other variables forms the argument arg to the MATLAB 
Bessel function. Since time varies for each time slice, arg varies for each time slice generating 
a filter per time slice. After generation each filter is stored to disk for use by AC-PROP.M. 
The variables base, NO, slices, and Step are also stored to disk for use by AC-PROP.M. 
The interested reader can find a detailed explanation and the source code in Appendix A. 
Graphical examples of the time- varying filters are include as Appendix B. 
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0 Spotiol frequency f x (cycle/m) 200 

Figure 3.2: Construction of rho shown graphically. 



3.2 Acoustic Propagation Module 

The script file, AC-PROP.M (ACoustic PROPagation), allows the user to specify the type 
and dimensions of the input excitation. It then takes the user’s chosen input function 
and simulates propagation as a function of time. AC -PROP computes the spatial impulse 
response p( 2 , j/, z, *), 

p{x,y,z,t) = TJ X \ jv {s(fz,fy)Hfx,fy,Z,t)} , (3.4) 

using the time-varying Bessel filters, produced by ACJFIL for h(f z , f y ,z,t). In Eq. 3.4, 
Tf xh {•} is the inverse spatial Fourier transform operator and the denotes a transformed 
function. 

To compute p(x,y,z,t), AC.PR0P first loads the variables passed from AC.FIL and 
then queries the user to make an input function choice. The user is given four input 
functions from which to choose: Circle, Table, Gaussian, and Bessel. The Table and Circle 
are equal-amplitude sources having the shape of a square and a circle, respectively (known 
in acoustics as the square and circular piston). The Gaussian and Bessel inputs are circularly 
truncated functions that have the indicated amplitude distribution within the circle (i.e., 
the amplitudes are spatially varying across the face of the source). Pursuant to the input 
function choice, the user is asked to input the diameter d of the truncating circle (or width of 
the truncating square w in the case of the Table function). Once the program was verified, 
a diameter d of 51 was used; this value is available as the default in the menus. In the 
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case of the Gaussian or Bessel input, the user is further requested to input the standard 
deviation sigma or a scaling factor a, respectively. These two function defining parameters 
were varied on a case-by-case basis. 

The variable shft-input holds the computed input array that was generated by the func- 
tion written to model the chosen excitation function (Appendix D). (The reader is reminded 
that shft prefix in a variable name indicates the array has the centered geometry as dis- 
cussed in Chapter II.) From shft-input , F -input is created by shifting (fftshift) shft-input 
to a comer geometry and then taking the two-dimensional spatial Fourier transform (fft2). 
As explained in Chapter II, the fftshift operation is necessary before the Fourier trans- 
form operation to obtain the correct phase relationship in the transform operation. With 
a shift back to the center geometry, the angular spectrum of the source 5(/ x ,/ y ), called 
shft-F-input in the program, is created for user viewing, if desired. The Bessel-function 
propagation transfer function h(/ x , f y ,z,t) from Eq. 3.1 must now be loaded so that the 
angular spectrum-propagation transfer function product sh on the right side of Eq. 3.4 
can be formed. The loading and multiplication process is repetitive since shft-F-input must 
form a product with the filter (or appropriate propagation transfer function) from each time 
slice. This repetitive multiplication is accomplished with a loop. 

The product of s and h (called shft -F .output). To find the desired result, the two- 
dimensional inverse spatial Fourier transform (ifft 2) must be taken. Before this can be 
done, F-output is formed by shifting shft-F.o utput from the centered geometry to the corner 
geometry. Executing the inverse transform of the product yields the output ( output ) which 
is then shifted to give shft.output for viewing. The array shft.output represents the output at 
the time slice that the loop is currently computing; shft.output does not depict the acoustic 
potential (or propagation pattern) through time; it only depicts the acoustic potential at a 
specific time. 

To produce a time history of the desired output, the center row (row NO) of shft.output 
is taken and placed in the array output.plot as the m-th column (m is the loop counter 
which relates directly to the time slice number, i.e., when m = 4, the computations for the 
fourth time slice are performed). This results in an array whose size is base x slices when 
Step zero-valued rows are added (that precede time slice NO row. The output examples 
are graphical interpretations of output-plot. Results generated in this manner are given 
in Chapter IV for all of the excitation functions. A detailed explanation of AC-PROP is 
provided in Appendix C. 

3.3 Program Summary 

The previous two sections gave an overview of the two program modules, AC.FIL and 
AC.PROP, including the code variable names and values assigned. What follows here is a 
summary of the steps that the program accomplishes. Step one is accomplished by AC.FIL. 
In this step the slices-Step filters to be used by AC-PROP are generated and saved. 

AC-PROP generates the user-specified excitation function s(i,y). Then the angular 
spectrum of the source s(f x ,f y ) is computed by taking the 2-D spatial Fourier transform 



16 



of s(x,y ). The product sh is computed for each time slice via a loop. In the loop, the 
inverse 2-D transform of the product forms an output for the specific time slice. The NO-lh 
(center) row of that output is then placed in successive columns of a new array to form the 
final output, p(x,0, 10, t). 

In the following chapter, examples of the final outputs are given. The excitations used 
for verification, as previously related, are the Table and Circle excitation functions. New 
values, as explain in this chapter, were then used for the simulation of the Gaussian and 
Bessel excitations. 
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Chapter 4 

Numerical Simulations 



In the previous chapter, a functional explanation of the two program modules was given 
including values assigned. The first section of this chapter reiterates the defining parameters, 
gives a brief explanation of each, and gives the parameter’s units. In the following section, 
the defining parameters are given the values used to verify correct operation of the program 
for the circular and square piston sources. The last section presents results for the non- 
piston Gaussian and the Bessel excitation functions. 



4.1 Defining Parameters 

A defining parameter is a parameter that delineates an aspect of the basic setup upon 
which all the remaining parameters or variables depend. In the work of this thesis there 
are two sets of defining parameters — those for the filter generation and those for generation 
of the excitation functions. The filter parameters are found at the beginning of AC.FIL. 
AC-PROP reads these parameters from a data file that is stored with the results of the 
filter calculations. 

The defining parameters found in ACJPIL include base, slices, Step, c, z, time.max, and 
rho.max. The first parameter base sets the dimensions of the base array, giving the number 
of sample points. The dimensions of the base array are, therefore, base x base where base is 
required to be a power of two (typically, 128). Making base a power of two allows MATLAB 
to use a high-speed radix-2 fast Fourier transform algorithm [16] to compute the spatial 
transforms. The next parameter slices is the number of time samples. Of these slices time 
slices, slices-Step slices require filters to be computed; the parameter Step is the number of 
leading-zero rows in the base x slices output array; as explained in the preceding chapter, 
this simulates the Heaviside step function. The parameters base, slices, and Step are unitless 
and are stored by AC _FIL in a file for use by AC-PROP. 

The remaining defining parameters of AC-FIL have units and are used only in the 
computations of the filters. The acoustic velocity in the medium, free-space in this case, 
is denoted by the parameter c having the units of m-s -1 . The source-to-reception point 
distance has the designation z with units of meters. The maximum time of propagation 
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Parameter 


Definition (units) 


base * 


Size of square base array 


slices m 


Number of time slices 


Step' 


Number of leading zero-rows 


c* 


Acoustic velocity in media (m/s) 


z* 


Distance, source-to-receiver (m) 


time-max * 


Maximum time of propagation (s) 


rho-max* 


Spatial radius of the filters (length -1 ) 


* passed to AC-PROP.M 



Table 4.1: Summary of the defining parameters in AC-FIL. 



time-max has units of seconds. The spatial-frequency radius of the filters rho.max (or rho ) 
has units of inverse length (i.e., m -1 , cm -1 , etc.). The unit of length depends on the area to 
be represented by the base array. These four parameters relate directly to Eq. 3.4 and are 
the parameters that dictate the diffraction properties of the filters. Table 4.1 summarizes 
the defining parameters found in AC.FIL. 

Another important set of defining parameters is the set that defines the user chosen in- 
put. These input defining parameters are entered by the user when requested by AC.PROP. 
Once the input function is chosen, the diameter of the truncation circle d is input. (The 
width of the table w (vice d) is input in the case of the Table excitation.) The parameters 
d (or w) are expressed as the number of points, out of base total points, that define the 
diameter (or width) of the function 

To transition from a diameter in terms of a number of points to an actual metric value, 
two equations were needed. The equations are 



and 



Ax = 



1 


(4.1) 


2 Pm ax 


k Ax 


(4.2) 


k 

2pmax 


(4.3) 



where Ax is the length of a segment, is the maximum spatial radius, k is the number 
of segments, and d is the diameter (or width w for the Table function). To determine the 
metric diameter, Eq. 4.1 was used to solve for Ax by setting /w* to 200 m -1 . This value 
resulted in a Ax of 2.5 x 10 -3 m or 2.5 mm. 

If the Gaussian excitation is selected, the user enters the value of the standard deviation 
<r, upon request. In a Bessel excitation selection the user enters the scaling factor a when 
requested. Table 4.2 gives a summary of the defining parameters used in AC_PROP. 
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Parameter 


Definition (units) 


base * 


Size of square base array 


slices * 


Number of time slices 


Step * 


Number of leading zero-rows 


c m 


Acoustic velocity in media (m/s) 


2* 


Distance, source-to-receiver (m) 


time-max* 


Maximum time of propagation (s) 


rho-max* 


Maximum spatial frequency (m -1 ) 


d 


Diameter of excitation function (sample points) 


w 


Width of Table excitation function (sample points) 


a 


Gaussian standard deviation . 


a 


Bessel scaling factor 


* passed from AC.FIL.M 



Table 4.2: Summary of the defining parameters used in ACLPROP. 



4.2 Program Verification 

In verifying the program output, two excitation functions were used, the Table and the Circle 
functions. The outputs generated by the program from these excitations were compared 
to the results found with the previous FORTRAN implementations [4, 3, 7] for validation. 
After a general explanation about the generation, formatting, and titling of the outputs, 
the two excitation functions are presented. The Table, the first verification function, is then 
discussed and a table of defining parameters is given. The second verification function, the 
Circle, follows with a similar discussion and table of defining parameters. 

4.2.1 Format of Results 

The graphical outputs for the two excitation functions used for verification, the Table and 
the Circle, were generated, formatted, and titled in the same manner. The outputs are for 
a source-to-receiver distance of z = 10 cm. Each output was for a given time slice and 
consisted of a 128x128 array of values. There were 64 time slices including 3 leading all-zero 
arrays which simulate the step function at the arrival of the wave ( t = z/c). The MATLAB 
programs contain optional commands to plot the output data using MATLAB’s plotting 
routines or to store the data in ASCII format for importation into the AXUM plotting 
program. (The AXUM program is more flexible than the MATLAB plotting routines.) 

Additionally, the 64 data arrays were combined with the Spyglass Data Utility into a 
128x128x64 array. Spyglass Transform and Dicer could be used to interactively pick the 
data slice of interest and to plot the results. 
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NAME 


VALUE 


DEFINITION 


base 


128 


Size of square base array 


slices 


64 


Number of time slices 


Step 


3 


Number of leading zero-rows 


c 


1500 


Acoustic velocity in media (m/s) 


z 


0.1 


Distance, source-to-receiver (m) 


time-max 


3.75 x 10 -4 


Maximum time of propagation (s) 


rho-max 


200 


Spatial radius of the filters (length -1 ) 


w 


23, 31 


Width of Table (samples) 



Table 4.3: Values of defining parameters for the Table input function used for program 
verification. 



4.2.2 Table Impulse Excitation 

The first excitation function to be run by the program was the Table function (for a square 
piston source). There were several reasons to use the Table as the first input; the Table 
function is an easy function to implement and the results could be readily compared to 
results found in the literature [3, 5, 6, 7]. Table 4.3 provides list of the defining parameters 
used. 

The values of base , slices , z, time-max, rho.max, and w chosen in Table 4.3 parallel 
those found in the literature used for validation. The acoustic velocity c of 1500 m/s is 
the velocity in water. To simulate the step function, a number Step of leading zero-rows 
was incorporated and arbitrarily set to three. The results were comparable to those in the 
literature and are presented (with the input functions) in Figs. 4.1 and 4.2. 

Here the widths w = 23 and w = 31 samples translate into metric widths of w = 5.5 cm 
and w = 7.5 cm, respectively. In the w = 23 case, k = 22 was the input in Eq. 4.3 and, in 
the w = 31 case, the input was k = 30 samples. Note that the x and y axes of the input 
functions range from 1 to 128, delineating a 128x128 base array. 

The outputs present the magnitude of the acoustic potential at an observation point 10 
cm from the source, p(x, 0, 10, <). A diffraction duration from the initial impulse ( t = 0) to 
t =time.max (150 ps) is represented as a function of radial distance; i.e., p(x, 0,10,0) to 
p(x,0, 10, 150) is represented. This gives the 3-D view of the general diffraction through 
time. The output images show several interesting features. 

The development of “tails,” explained in terms of edge waves [3, 18], can be seen in 
the 3-D images. Also of note are the overshoots, having a maximum magnitude of 2.11 
(both cases), and the undershoots (difficult to see in these two cases); these are due to the 
additive nature of the interference patterns of the waves originating on the edges of the 
discontinuous source. 
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(o) 




(b) 



Figure 4.1: Table spatial input and time-space output for w = 23 samples (w = 5.5 cm) at 
z — 10 cm. 
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Figure 4.2: Table spatial input and time-space output for w = 31 samples (w = 7.5 cm) at 
z — 10 cm. 
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Figure 4.3: Table output for d = 31 samples. Dicer representation. 



3D outputs for Table excitation 

The Spyglass Dicer program allows selected three-dimensioned representations of the propa- 
gation output. Figure 4.3 shows a rectangular cube set into one-fourth of the data volume. 
The x and y axes represent the spatial variables and range from 0 to 127 samples. The time 
axis represents the 64 time slices of data that were computed. To elongate the time axis, 
the Dicer program was used to increase the number of slices by a factor of 4x (to a total of 
256 slices); linear interpolation is used to compute the values of the expanded slices lying 
between the original slices. Due to the nature of the interpolation used, only the expanded 
slices ranging from 0 to 252 are shown on the time axis. The three visible surfaces of the 
rectangular cube show the calculated data in three planes. The three planes are located at 
x = 64, y = 65, and time = 252. An additional plane is located at time = 16 to show the 
shape and size of the source (since the theory predicts that the field at time = 16 expanded 
samples (when t = z/c] duplicates the source excitation). 

An alternative Dicer representation of the field is shown in Figure 4.4. Here, a horizontal 
slice of the data is placed at y = 65 and five vertical slices of the data are located at x = 



24 





- 0.1455 



&ne 



Amplitude 



240 



2.070 



127 



y 



o 






Figure 4.4: Table output for <i = 31 samples. Alternative Dicer representation. 
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16, 64, 132, 200, and 252. (The values were arbitrarily chosen, but were restricted to be a 
multiple of four to avoid showing interpolated values of data.) 

Appendix E contains Dicer representations for the other Table and Circle excitations 
that were studied. 

4.2.3 Circle Impulse Excitation 

The defining parameters for the Circle excitation are the same as those introduced in Ta- 
ble 4.3 with the exception that only the d = 31 case is presented. A diameter d = 31 samples 
translates into a metric diameter of d = 7.5 cm. Again the results are comparable to those 
found in the literature [5, 7], Figure 4.5 gives the input function and the ensuing outputs. 
As in the Table case, the base array is a 128x128 array with the propagation pattern formed 
by successive p(x,0, 10,2) time slices. The results in Fig. 4.5 for the Circle excitation are 
much the same as those for the Table in Fig. 4.2. The “tails,” however, shown in the 3-D 
image of Fig. 4.5 are rounded instead of cornered as in the Table output. Though the Table 
output holds its magnitude for a longer time, the maximum for the Circle output is greater 
at 2.21. (This value was read from the array; it is difficult to obtain quantitative informa- 
tion from the 3-d plots.) Also the drop off from the maximum is steeper for the Circle . The 
greater maximum and steeper drop off are due to the equal distance of all edge points from 
the center. This same geometric influence also accounts for the Circle holding the input 
value for a shorter duration. Just as the interference patterns added to a maximum greater 
than the input, the interference patterns also combine to give a more negative minimum. 
The negative undershoot was present for the Table ; however, it has a magnitude five times 
greater for the Circle . 

4.3 Other Input Excitations 

Having checked the performance of the technique with piston sources, other sources with 
nonuniform spatial excitations were investigated The first is a circularly truncated Gaussian 
distribution function. Following the Gaussian, a circularly truncated Bessel profiled function 
is examined. These two excitation function outputs are generated and formatted the same. 

4.3.1 Gaussian Distributed Excitation 

Though the Gaussian has been investigated before [3]- [7], it has not been studied as 
a 128x128 array. The defining parameters for the case studied are listed in Table 4.4. 
Figure 4.6 shows the input and resulting output. 

The Gaussian excitation function has been normalized by the maximum value of the 
computed Gaussian (see the Gaussian source code titled CRCGAUS.M in Appendix D). 
This normalization is shown in the input image of Fig. 4.6 by the maximum amplitude of 
one. Displayed as a base x base array, this input image has a standard deviation of a = 5 
and a 1/e point of 10.17 samples from the center (metric equivalent of r = 2.54 cm). The 
diffraction of this input is shown in Fig. 4.6. 
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Figure 4.5: Circle input excitation and output for d — 31 samples. 
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NAME 


VALUE 


SUMMARY 


base 


128 


Size of square base array 


slices 


64 


Number of time slices 


Step 


3 


Number of leading zero-rows 


c 


1500 


Acoustic velocity in media (m/s) 


z 


0.1 


Distance, source-to-receiver (m) 


time-max 


375 x 10 -6 


Maximum time of propagation (s) 


rho-max 


200 


Spatial radius of filters (length -1 ) 


d 


51 


Diameter of excitation function (samples) 


a 


5 


Gaussian standard deviation 



Table 4.4: Defining parameters for Gaussian excitation case. 



NAME 


VALUE 


SUMMARY 


base 


128 


Size of square base array 


slices 


64 


Number of time slices 


Step 


3 


Number of leading zero-rows 


c 


1500 


Acoustic velocity in media (m/s) 


z 


0.1 


Distance, source-to-receiver (m) 


time-max 


375 x 10 -6 


Maximum time of propagation (s) 


rho-max 


200 


Spatial radius of filters (length -1 ) 


d 


51 


Diameter of excitation function (samples) 


a 


0.25 


Bessel scaling factor 



Table 4.5: Defining parameters for Bessel excitation case. 



The 3-D image shows a diffraction pattern that is well established by time t = time. max , 
forming two spreading “tails.’ 1 The “tails,” as well as the rest of the Fig. 4.6. diffraction 
pattern, are smoothly rounded. This rounding is the result of the continuity of the Gaussian 
distribution. A discontinuity, as in the previous two excitation shapes, results in a char- 
acteristic over and undershoot of the maximum and minimum inputs. Again, the results 
for this Gaussian excitation conformed to those found in the literature [3, 4, 5, 6, 7]. The 
Bessel excitation was then run and the results compared to those of the Gaussian. 

4.3.2 Bessel Excitation 

Results for the Bessel excitation produced by CRCBES.M were generated, as previously 
discussed, for the set of defining parameters listed in Table 4.5. 

The resulting input and output are shown in Fig. 4.7. Since the output is formed by 
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(a) 




Figure 4.7: Bessel-profile input and output for a = 0.25. 
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taking successive p(x,0,10,t) vectors, three peaks appear in the 3-D image. These three 
peaks correspond to the center peak and the points on the two crests directly adjacent 
to the center. As a result of having three peaks, three “tails” are present. The “tails,” 
however, are smooth because a Bessel function is a continuous function. Also worth noting 
in the 3-D image is the simulated step function, more visible here due to the oscillations of 
a Bessel. 

Comparing the Gaussian and Bessel outputs, it is seen that the Gaussian’s magnitude 
retention is slightly longer than that of the Bessel. This is due to the more gradual decrease 
in magnitude vice the steeper decrease required of the Bessel so that it can become negative. 
Still no hard conclusions can be drawn without further analysis. 
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Chapter 5 

Summary 



This report presented a MATLAB implementation of a Fourier approach to acoustic wave 
propagation. A mathematical development using linear systems that found the acoustic 
potential from an arbitrary spatial and temporal source was reviewed. In the mathematical 
development, it was shown that the Green’s function solving the appropriate wave equation 
and satisfying the boundary conditions is the total impulse response of the system. Through 
double and triple convolutions, the acoustic potential could be found for any source sepa- 
rable in time and space. Use of the 2-D spatial Fourier transform, however, translated the 
convolution to multiplication in the spatial frequency domain. This made the MATLAB 
implementation easier. 

After an overview of MATLAB and the graphics program AXUM, a functional descrip- 
tion of the program was furnished. The program modules AC.FIL and AC.PROP both 
made use of symmetry. AC-FIL generated the time-varying filters, the most time con- 
suming process, while AC-PROP accomplished the remaining computations making use of 
MATLAB ’s “fft2” function. Details of both modules as well as the source code have been 
included in the Appendix A. 

Several examples were delineated in the body of this thesis. First, the Table and the 
Circle were presented as the program verification excitations; the results conformed to those 
found in the literature. Then two newer excitation functions, the Gaussian and the Bessel, 
were presented. 

The underlying result was an accurate and efficient computer implementation of the 
linear systems approach to ultrasonic wave propagation. The efficiency was derived from the 
modularization of the program so that consecutive runs could be made without recomputing 
the most time consuming portion, the filters. Also, the use of MATLAB ’s “fft2” function 
bypassed tedious and time-inefficient convolution integrals. Finally, both modules made 
use of symmetry by computing only one quadrant of data which was then manipulated into 
the remaining quadrants. An advantage to using MATLAB was the ease of expansion that 
could be accomplished with the program. 

The work of this report concentrated on a source with rigidly baffled boundary condi- 
tions and a lossless media. Cases that include free space and resiliently baffled boundary 
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conditions as well as lossy media, linear and quadratic lossy media, could be incorporated. 
A few facts worth noting here are that the free space and resilient baffle boundary conditions 
can be expressed in terms of the rigid baffle case [4] and that the lossy media and lossless 
media transfer functions are interdependent [6]. Furthermore, new excitation functions, 
such as a phased array or a focused source [2], could also be incorporated. Improvements 
are needed in the area of analysis such as the Gaussian versus Bessel propagation compar- 
ison and extending the technique to sources that are not time and space separable such as 
new non-diffracting waves [19]. 
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Appendix A 

Source Code for AC_FIL.M 



The following is the source code used to generate the time-varying filters as discussed in 
Chapters II and III. AC-FIL.M was written in block format with each block headed by 
a descriptive comment to explain the block’s function. The code includes many optional 
instructions indicated by the leading %* symbol. Deleting this symbol will enact that 
line of code on succeeding program runs which, in turn, varies the output. The outputs 
necessary to a successful run of AC_PROP.M are the files, ACbasex(m + Step).MAT (where 
m is an index number from 1 to 61). and the file, AC.FIL.MAT. For example, the file, 
AC128x08.MAT, contains the data for the propagation spatial filter in a 128x128 array for 
the fifth time slice (since Step = 3). The file, AC-FIL.MAT, contains the parameters needed 
in AC-PROP.M. 

Code is provided for both the DOS version and the Sun workstation version of MAT- 
LAB4. There are two primary differences in the versions. First, the paths to disk storage 
are different, to reflect the path setup on the two host machines. Secondly, the commands 
to obtain a hard-copy version of MATLAB4’s graphics differed. In the DOS version, the 
user prints copies from the menu associated with the graphics window or with a print com- 
mand. In the SUN version, a hard copy can be printed only by a print command. Each 
method is included in the file text; the user selects the appropriate command by removing 
the “the desired lines. 

ACJFIL.M SOURCE CODE 



7, **** AC.FIL.M **** 

7.7. 

7.7. This program generates an Acoustic Propagation Transfer 

7.7. Function, a time varying spatial filter, for use in 
7.7 AC.PROP.M to simulate acoustic wave diffraction. 

7.7. William H. Reid December 1992 

7.7. Modified for MATLAB4 and Sun. 9/93 JPP 

clear; 7. Clear MATLAB 
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tml=clock; ; 



'/, Start timer clock 



rho.max = 200; 



c = 1500; 
z = 0.1; 

time.max = 3.75e-4; 



base - 128; 

NO = (base/2)+l ; 
slices = 64; 

Step = 3; 



7, Size of square base array. 

'/. Defines center of square base array. 

'/, Number of time slices. 

'/. Number of leading zero time slices, 
simulates the step function. 

'/, Velocity of the acoustic wave, (m/s). 

'/, Distance to the observation plane, (m) . 
'/. Maximum time of propagation, 

'/ time at the final time slice, (sec). 

'/ Spatial radius of the filter, 

7, [sqrt (rhox~2 + rhoy~2] . (per length) 



•/.y. Initialize arrays to reduce processing time, 
shft.filter = zeros(base); temp = zeros(NO); 
arg = zeros(NO); rhom = zeros (NO, 1); 
rho = zeros(NO); time = zeros (slices-Step , 1) ; 

VL Generate "slices-Step" time slices between z/c and time.max. 
time = linspace(z/c, time.max, slices-Step) ; 

y.y, Choose directory to store data. 

cd i : /ac_ prop/filters '/, PC version 

'/,* cd /home2/powers/M„f iles/ac„prop/f ilters */, SUN version 
y.y. Save variables necessary for passing to AC_PR0P.M in AC_FIL.MAT. 
save ac.fil base NO slices Step c z time.max rho.max; 

y//. Generate N0-1 values of "rhom 11 from 0 to rho.max. 
rhom = linspace(0,rho.max ,N0-1) ; 

y//. Add additional increment to rhom to compensate for off-center 
y.y, orientation of the final base x base matrix, 
rhom = [rhom (rhom(N0-l)+rhom(2) )] ; 

y//. Create two NO x NO arrays of rho values for function evaluation. 
[rhox,rhoy] = meshgrid (rhom, rhom) ; 

y//. Calculate “rho", an NO x NO matrix of radial distances for use in 
•/.•/. the argument to the Bessel function within the loop. 
rho= sqrt(rhox.~2 ♦ rhoy.~2); 
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mmmm startloop murnu 

Generate "slices-Step" filter arrays, the filter at each time slice, 
for m = 1: (slices-Step) 

fprintf (’7.3. Of * ,m) ; */, Display m value for user progress report. 



7.7. Create an NO x NO array of argument values for the Bessel function, 
arg = rho * sqrt(c"2*time(m)~2-z~2 ); 



Evaluate the zero order Bessel for each argument value; 
creates an NO x NO axray called "temp", 
temp = 2*besselj (0,arg) ; 



Create shft_filter matrix containing the spatial filter. 

7.7. (Done by flipping "temp" into all quadrants.) 
shft_filter(NO:base,NO:base) * temp(l :N0-1 ,1 :N0-1) ; 
shft_f ilter(NO:base, 1 : NO- 1 ) = f liplr (temp(l :N0-1 , 1 :N0-1) ) ; 
shft.f ilter(l :N0-1 , 1 :N0-1) = temp(NO:-l :2,N0:-1 :2) ; 
shf t_f ilter(l :N0-1 ,N0:base) = flipud(temp(2:N0,2:N0)) ; 



Save shft.filter in a file named "a(base)x(m+Step) .mat" , 
cd i : /ac_prop/f ilters '/, PC version 

'/.* cd /home2/powers/M_files/ac_prop/f ilters 7. SUN version 



if (m+Step) < 10, 

7.7. MATLAB format 

evaKC’save a’ ,int2str(base) , ’x0’ , int2str (m+Step) , ’ shft.filter m’ ] ); 

7.7. ASCII format 

7.* eval([’save a’ ,int2str(base) , ’x0’ ,int2str (m+Step) , ’ .dat shft.filter. 
’/,* /ascii’]); 

else 

7.7. MATLAB format 

eval([’save a’ ,int2str(base) , ’x’ ,int2str (m+Step) , ’ shft.filter m’ ] ); 

7.7. ASCII format 

’/,* eval([’save a’ ,int2str(base) , ’x’ ,int2str(m+Step) , ’ .dat shft.filter . . 
'/,* /ascii’]); 

end 



end 

7.7.7.7.7.7.END LOOP 7.7.7.7.7.7.7.77. 

toe; 7. Stop timer clock 

time.elapsed « etime(tm2 ,tml)/60 7. Compute & display execution time 
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Appendix B 

Examples of the Time-varying 
Propagation Bessel Filters 



This appendix contains examples of the time-varying filters (Eq. 3.1 on page 12) produced 
by AC-FIL for the parameters described in the text. Only a few of the 61 total slices are 
shown. The first slice, at t = z/c, is a plane with amplitude two; this is because the input 
argument to the Bessel function is zero giving an output value of one. This uniform plane 
produces an output for time slice one that is a scaled replica of the input. The remainder of 
the filters illustrate the time variance and how the filters collapse inward with time. Each 
filter is a 128x128 array representation in the spatial frequency domain. 
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(a) 





(c) 



Figure B.l: Propagation filter at time slices 1, 2, and 5. 
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(a) 






Figure B.2: Propagation filter at time slices 10, 20, and 30. 
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Appendix C 

Source Code for AC_PROP.M 



The following is the source code for AC-PROP.M that was used to produce various outputs, 
including those discussed in Chapter IV. AC-PROP was written in blocks with each block 
headed with a descriptive comment to explain the block’s function. This also makes it 
easy to follow the computation from inputs to output. Inputs to AC-PROP are imported 
from the file AC.FIL.MAT and the files A basexm + Step. MAT (m is an index number 
ranging from 01 to 61). The program solicits user input to determine the input excitation. 
AC-PROP then computes the output. 

The format of the output can be changed by the user. Before running the program, the 
user can remove the optional comment markers indicated with %% to produce and/or view 
the output in the desired format. This allows the data to be saved and exported in ASCII 
format for use with other graphics programs (such as AXUM, see Chapter II). 

AC.PROP.M SOURCE CODE 



r/.m **** AC.PROP.M **** 

n 

This script file performs transient -wave acoustic propagation 
*/,*/, simulations. It uses the time-varying spatial filters called 
*/,'/, "shft.f ilter" found in the filters directory under "a(base)x(m) .mat" 
*/,'/, to compute the acoustic propagation fields. The "shft.f ilter" 

'/,*/, data are generated by AC.FIL.M. 

*/.*/. 

William H. Reid, December 1992 
*/.*/. Modified for MATLAB4 and Sun. 9/93 JPP 

clear; clc; '/, Clear MATLAB 

tic; */, Start timer 

format compact '/, Set compact format for screen display. 
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VI Load the parameters generated by last run of AC.FIL.M. 

VI PC and SUN directories have different names 

cd h:\ac_prop\filters */, PC version 

'/,* cd /home2/powers/M.f iles/ac. prop/filters '/, SUN version 
load AC.FIL 

*/,'/, Display the size of the square base array. 
dispCbase is the width of the array. 1 ); dispC 1 ) ; base 

H Generate the INPUT function from user interface. 

input. func = menu( 'Choose the shape of input function. ', 'Circle' , ... 

'Square' , 'Truncated Gaussian' , 'Truncated Bessel'); 
if isempty (input. func) ; input. func = 3; end '/, Default to Gaussian input, 
if input. func == 1, */, Cicle input 

name=' c' ; 

dispCYou have chosen a truncated circle input.') 
d = input ( 'Please enter an ODD diameter, [51], d = '); 

if isempty(d); d - 51; end */, Default diameter: 51 samples 

shft. input = crcle(d,base) ; '/. Create Circle input 

'/, Name a file to hold the input function. 

InFile.Name = [name, 'i' , int2str (base) , 'x' ,int2str(d)] ; 

'/, Name a file to hold the output. 

OutFile.Name = [name, 'o' , int2str(base) , 'x' ,int2str(d)] ; 
elseif input. func == 2, '/, Square input 

name= ' s ' ; 

dispCYou have chosen a truncated square input.') 
w = input ( 'Please enter an ODD width, [11], w = '); 

if isempty(w); w = 11; end '/, Default width: 11 samples 

shft. input = table(w,base) ; '/, Create input 

d - w; 

*/ Name a file to hold the input function. 

InFile.Name = [name, 'i' , int2str (base) , 'x' ,int2str(d)] ; 

'/, Name a file to hold the output. 

OutFile.Name = [name, 'o' , int2str (base) , 'x' ,int2str(d)] ; 
elseif input .func 3, '/, Gaussian input 

name='g' ; 

dispCYou have chosen a truncated Gaussian input.') 

sigma = input ('Please enter the standard deviation, [10], sigma = '); 

if isempty (sigma) ; sigma = 10; end '/ Default sigma: 10 samples 
d = input ( 'Please enter an ODD diameter, [51], d = '); 

if isempty(d); d s 51; end '/, Default diameter: 51 samples 
shft. input = crcgaus(sigma,d,base) ; '/. Calculate input 
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'/, Name a file to hold the input function. 

InFile.Name = [name, ’i’ , int2str(base) , ’x’ ,int2str(d)3 ; 

7. Name a file to hold the output. 

OutFile.Name = [name, ’o’ , int2str(base) , ’x’ ,int2str(d)] ; 
elseif input.func == 4, 7. Bessel input 

name='b’ 

disp(*You have chosen a truncated Bessel input.’) 
a = input( ’Please enter a width scaling f actor , [0 . 3125] , a = ’); 

if isempty(a); a = 0.3125; end 7. Default: 0.3125 
d = input( ’Please enter the ODD diameter, [51], d = ’); 

if isempty(d); d = 51; end 7. Default: 51 samples 

disp(’Please wait. This calculation takes a while ’) 

shft.input = crcbess(a,d,base) ; 7. Calculate input 

q = a * le4; 

7. Name a file to hold the input function. 

InFile.Name = [name, ’i’ , int2str(base) , ’x* ,int2str(q)] ; 

7, Name a file to hold the output. 

OutFile.Name = [name, ’o’ ,int2str(base) , ’x’ ,int2str(q)] ; 

lse 

disp(’ *); 

disp( ’Incorrect Excitation Function Selection!’); 
error ( ’Restart AC_PR0P...to try again.’); 

end 

7.7. Save the input function placed in "InFile.Name" in ascii format 

7.7. for use with other graphics software. 
eval([’save ’ ,InFile_Name, ’ .dat shft.input /ascii’]); 
clear InFile.Name; 7. Remove "InFile.Name" from memory. 

7.7. Compute the sample spacing in x and y (delta.x) and the spacing in 

7.7. time (delta.t) . 

delta.x = (base-2)/ (2*rho_max) ; delta.t = time_max/(slices-Step) ; 

7.7. Create X and T vectors (1 x base). 

X = linspace(-(NO-l)*delta_x, C(base-NO)-l)*delta_x,base) ; 

T = linspace(-3*delta_t, time.max, slices) ; 

7.7. Display the input function. 

7.* mesh (X , X , abs (shf t.input) ) ; t itle ( ’ | SHFT.INPUT | ’ ) ; 

7.* axis( [ X(l) X(base) X(l) X(base) ... 

7.* min( [ 0 min(min(abs(shft. input))) ] ) max(max(abs(shft_input))) ]) 

7.* xlabelC’x (cm)’); ylabel(’y (cm)’); zlabel( ’amplitude’) ; view(52.5,30) ; 
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*/,'/, SUN: Save input figure as eps file. 

7* cd /home2/pouers/M_files/ac_prop/data 

77 PC: save figure from figure window. 

7* cd h:/ac_prop/ 



Save input plot as EPSF file with bitmap preview image. 

’/,* disp(’ ’); disp(’Saving input plot as EPS file ’); disp(’ 

7* eval( [’print -deps -epsi ’ .name, ’.in.b’ ,int2str (base) , ... 

7* ’_d’ ,int2str(d)] ) ; 

7 * disp(’ ’); disp(’Saving input plot as EPS file ’ ) ; disp(’ 

7* eval([’print -deps -epsi ’ .name, ’.in.b’ ,int2str (base) , ... 

'/.* ’ _d’ , int2str(d)] ) ; 



’): 

’); 



7 . 7 . Shift "shft.input" from centered geometry to corner geometry and take 

7 . 7 . 2-D FFT to produce F.INPUT. 

F.input = fft2( f ftshif t(shft.input) ); 

clear shft.input; 7 . Free RAM. 

7 . 7 . Shift F.input in preparation of multiplication with PROP.m. 
shft.F.input = ff tshif t(F.input) ; 
clear F.input; 7 . Free RAM. 

7 . 7 . Element by element arrray multiplication of the transfer function 

7 . 7 . filter in "PROP.m" and "Fshft.input . " 
disp( ’Performing array multiplication....’); 



7 .* cd /home2/powers/M_f iles/ac.prop/data 7 . SUN version 
cd h:/ac_prop/data 7 . PC version 



7 . 7 . Save "Step" arrays full of zeros for first array values. 

7 . 7 . 7 . 7 . 7 . 7 . 7 . 7 . 7 . 77 . 7 . 7 . START FIRST LOOP 77 . 7 . 7 . 7 . 7 . 77 . 7 . 7 . 7 . 7 . 7 . 7 . 7 . 



for m=l :Step 

shft.output = zeros(base); 

7 . 7 . Save the (base)x(base) array "shft.output" in an ASCII file whose 

7 . 7 . name is " (input.func) .OUT. (m) .dat" . 

eval([’save ’ ,int2str(input_func) , ’_0UT_O’ ,int2str(m) , ... 

’.dat shft.output /ascii’]); 

clear shft.output; 7 Get ready for next pass. 

end 



77 . 77 . 7 . 77 . 7 . 77 . 7 . 7 . 7 . END FIRST LOOP 77 . 7 . 7 . 7 . 7 . 7 . 7 . 77 . 7 . 7 . 7 . 7 . 7 . 
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mxmmm start second loop 77.7.77.7.7.7.7.7.777.77 

for m = 1 : (slices-Step) 

fprintf (’7.2.0f , ’ ,m) 7. Display "m" value for progress report. 



7.7. Give the variable "filename!." the name of the file containing 
7,'/, "shf t_f ilter" and then load that file, 
if m < 10-Step, 

filenamel = [’a’ ,int2str(base) , ’xO' , int2str(m+Step)] ; 
else 

filenamel = [’a* , int2str(base) , ’x’ ,int2str(m+Step)] ; 
end 

7* cd /home2/powers/M_f iles/ac.prop/f ilters 7, SUN version 
cd h:/ac_prop/filters 7. PC version 

eval([’load '.filenamel] ); 



7.7. Multiply the filter by the shifted transform of the input. 
shft_F_output = (shft.filter .* (shft_F_input)) ; 
clear shft.filter; 7. Free RAM 



7.7. Shift "shft_F_output" from centerd geometry to corner geometry 

7.7. ("F.output") and take inverse transform to produce "output." 

output = ifft2( fftshift(shft_F_output) ) ; 

7.7. Shift "output" to center geometry ("shft.output") . 
shft.output = fftshift (output) ; 

7.7. "shft.output" is the centered geometry version of the diffracted 

7.7. wave at time slice "m" . 

clear shft.F.output output 7. Free RAM. 

7.7. Place the transposed "NO" (center) row of "shft.output" in the 

7.7. (m+Step)-th column of "output _plot . " 

7.7. This creates a (base) x (m+Step) array whose columns show a slice 

7.7. of the diffracted wave. 

output_plot(l:base, m+Step) = (shft_output(NO, 1 :base)) ’ ; 

7.7. Save the (base)x(base) array "shft.output" in an ASCII file whose name is 

7.7. "( input _f unc) _ OUT. (m+Step) .dat". 

7.7. A Gaussian, for example, would be G_0UT_12.dat on the 9-th loop 

7.7. interation (Step = 3) . This is optional for graphics use by 
7.7 other software. 

7* cd /home2/powers/M_f iles/ac.prop/data 7 SUN version 
cd h:/ac_prop/data 7. PC version 

if m < 10-Step, 
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eval([’save ’ , int2str( input. f line) , ’.OUT.O’ , int2str (m+Step) , ... 
’ .dat shft.output /ascii’]); 
else 

eval([’save ’ , int2str(input_func) , ’_0UT_’ ,int2str(m+Step) , ... 
’.dat shft.output /ascii’]); 
end 

clear shft.output; '/, Get ready for next pass, 

end 

mmmmnw. End loop xxxxxxxxxxxxxxx 



'/,*/, Save the contents of "output.plot" in a MATLAB file 
*/,'/, and as an ascii file (optional). 

'/,* cd /home2/povers/M_f iles/ac_prop/data '/, SUN version 
cd h:/ac_prop/data */, PC version 

eval([’save ’ .OutFile.Name, ’ output.plot’]) 
eval([’save ’ .OutFile.Name, ’ .dat output.plot /ascii’]); 



’/,*/. Display output view #1. 

7.7. figure; mesh(T,X , abs (output.plot) ) ; title (’ I OUTPUT | ’) ; 

7.7. axis( [ T(l) T(slices) X(l) X(base) ... 

*/,’/, min( [ 0 min (min(abs (output.plot))) ] ) ... 

7.7. max(max(abs(output_plot))) ]) 

7.7. xlabel(’time (s ) ’ ) ; ylabel(’y (cm)’); zlabel(’ abs (output) ’) ; 
*/.•/. view (52. 5,30); 



'/.'/, Save output figure as eps file. 

*/.* cd /home2/pouers/M_f iles/ac.prop/data '/. SUN version 
cd h:/ac_prop/data '/, PC version 

'/,* disp(’ ’); disp(’ Saving output figure as EPS file ’); disp(’ ’); 

'/,* eval([’print -deps -epsi ’ .name, ’.outl.b’ ,int2str(base) , ... 

'/,* ’_d’, int2str(d)] ) ; 



77. Create a different view of output 

'/,* figure; mesh(T,X,abs(output_plot)) ; title(’ I0UTPUTI ’) ; 

•/.* axis ( [ T(l) T(slices) X(l) X(base) ... 

'/,* min( [ 0 min(min(abs (output.plot))) ] ) ... 

'/,* max(max(abs(output_plot))) ]) 

'/,* xlabeK’time (s)’); ylabel(’y (cm)’); zlabel(’ abs (output) ’) ; 

'/,'/, Save second output figure as eps file. 

'/.* disp(’ ’); disp(’Saving figure as EPSF file ’); disp(’ ’); 

'/* eval([’print -deps -epsi ’ .name, ’_out2_b’ ,int2str(base) , ’_d’ , ... 
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int2str(d)] ) ; 



*/.♦ 

elapsed_time = toc/60 '/, Stop timer; display elapsed time. 
disp( 'minutes ’ ) 
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Appendix D 

Source Code for Input excitations 



The following source code is for the input excitation choices given to the user in AC.PROP. 
Each is written as a MATLAB function and can be used independently of AC.F1L or 
AC.PROP. The input excitations are the uniform square ( TABLE(w,base ) ), the uniform 
circle ( CRCLE(d,base ) ), the circularly truncated Gaussian ( CRCGAUS(sigma,d,base )), and 
the circularly truncated Bessel ( CRCBES(a,d,base )) where w is the width of the square, d 
is the diameter of the circle, a is the standard deviation of the Gaussian distribution, a is 
a scaling factor, and base is the size of the base array. 

TABLE.M SOURCE CODE 



function Y = table(w.base) 

V. TABLE.M: Y = table(w.base) 

'/.Program for generating a uniform square excitation function. 

'/, December 1992 William H. Reid 
'/. Based on TABLE.M by JG Upton 
7 . 

7, w is the WIDTH of the table. (ODD integer) 

7% base is the WIDTH of the square base. (EVEN integer) 

7, Example: z - table(33,64) ; 

'/ Check that w is an odd integer, 
if rem(u,2) < 0.1; 

error(’The width of the table must be an ODD integer.’); 

else; 

end; 

'/, Check that base is an even integer, 
if rem(base,2) "= 0.0; 

error(’The width of the square base must be an EVEN integer.’); 
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else; 

end; 



NO = (base/2) + l ; 
wO « ceil(w/2); 



'/, NO is the base array’s center. 

*/, vO is the mid — point of the table. 



Y = zeros(base) ; 
temp » zeros(NO-l); 



7. Initialize matrices to reduce 
7. processing time. 



temp(l:vO,l:vO) = ones(vO); 



7, Set amplitude to one. 



7. Generate the entire base x base input function array. 
Y(baseO :base ,baseO :base) = temp; 

Y(2:base0,base0:base) = rot90(temp); 

Y (2 :N0 ,2 :N0) = rot90(temp,2) ; 

Y(N0:base,2:N0) = rot90(temp , 3) ; 

'/, To test input distribution: mesh(Y) 



CRCLE.M SOURCE CODE 



function Y = crcle(d.base) 

7. CRCLE.M: Y = crcle(d.base) 

7, Program for generating uniform circular excitation functions 
7. December 1992 William H. Reid 
7. Based on CRCLE.M by JG Upton 



7. d is the DIAMETER of the circle. (ODD integer) 

7. base is the WIDTH of the square base. (EVEN integer) 

7. Example: z = crcle(33,64) ; 

7. Check that d is an odd integer, 
if rem(d,2) < 0.1; 

error (’The diameter of crcle function must be an ODD integer.’); 

else ; 

end; 

7. Check that base is an even integer, 
if rem(base,2) "= 0.0; 

error (’The width of the square base must be an EVEN integer.’); 



7. 
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else; 

end; 



NO * (base/2)+l; 
r = d/2; 



'/. NO is the base array’s center. 
'/. r is the circle’s radius. 



’/, Initialize matrices to reduce processing time. 

Y = zeros(base) ; 
temp = zeros(NO-l); 

'/. Set amplitude to one inside the circle’s radius, 
for m = 1 :r+l ; 

for n = 1 :r+l ; 

if sqrt((m-l)“2 + (n-l)“2) <= r; 

temp(m,n) = 1; 

end; 

end; 

end; 

*/, Generate the entire base x base input function. 
Y(baseO :base ,baseO :base) = temp; 
Y(2:base0,base0:base) = f lipud(temp) ; 

Y(2:N0,2:N0) = rot90(temp ,2) ; 

Y (NO:base ,2 : NO) = fliplr (temp) ; 

'/, To test input function distribution: mesh(Y) 

CRCGAUS.M SOURCE CODE 



function Y = crcgaus(sigma,d,base) 

*/, CRCGAUS.M: Y = crcgaus(sigma,d,base) 

*/, Program for generating circular Gaussian excitation functions. 
*/, December 1992 William H. Reid 
'/. Based on CRCGAUS.M by JG Upton 



'/, sigma is the STANDARD DEVIATION of the Gaussian function. 
*/. d is the DIAMETER of circle. (ODD integer) 

’/, base is the WIDTH of the square base. (EVEN integer) 

'/. Example: z = crcgaus(l2,33,64) ; 



•/. 



mu=0; 



'/.mu is the mean of the Gaussian function. 
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'/, Check that d is an odd integer, 
if rem(d,2) < 0.1; 

error(*The diameter of circle function must be an ODD integer.’); 

else; 

end; 

'/, Check that base is an even integer, 
if rem(base,2) "= 0.0; 

error(’The width of the square base must be an EVEN integer.’); 

else; 

end; 

NO = (base/2) + 1; '/, NO is center of the array, 

r * d/2; ’/, r is the radius of the truncating circle. 

'/, Initialize the matrices to reduce processing time. 

Y = zeros (base) ; 
temp * zeros(NO-l); 

*/. Compute the amplitude for the Gaussian distributed circle, 
for m * 1 : (d+l)/2; 

for n = 1 : (d+l)/2; 

x = sqrt((m-l)*2+(n-l)~2) ; 
if x <= r; 

temp(m,n) = (l/(sqrt(2*pi)*sigma))*exp(-((x-mu)~2)/ . . . 
(2*(sigma~2))) ; 

end; 

end; 

end; 

'/. Generate the entire base x base input array. 

Y(baseO:base,baseO:base) * temp; 

Y(2:base0,base0:base) = flipud(temp) ; 

Y(2:N0,2:N0) = rot90(temp ,2) ; 

Y(N0:base,2:N0) = fliplr(temp) ; 

Y = Y ./ (max(max(Y))) ; '/. Normalize the Gaussian distribution to one. 

*/, To test and view the input function: mesh(Y) 
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CRCBESS.M SOURCE CODE 



function Y = crcbess(a,d,base) 

'/. CRCBESS.M: Y = crcbess(a,d,base) 

'/ Program for generating circular Bessel excitation functions. 

'/, December 1992 William H. Reid 
'/, Based on CRCBESS.M by JG Upton 
V. 

I a is the WIDTH SCALING FACTOR. 

*/. d is the DIAMETER of the circle. (ODD integer) 

'/, base is the WIDTH of the square base. (EVEN integer) 

*/. Example: z = crcbess(l,33,64) ; 

*/. Check that d is an odd integer, 
if rem(d,2) < 0.1; 

error(’The diameter of the circle must be an ODD integer’); 
else; 
end; 

'/, Check that base is an even integer, 
if rem(base,2) “= 0.0; 

error(’The width of the square base must be an EVEN integer’); 
else; 
end; 

NO = (base/2) + l; '/ NO is the center of the array, 

r = d/2; '/. r is the radius of the circle. 

Y = zeros (base); */, Initialize the arrays to reduce 

temp = zeros(NO-l); */, processing time. 

'/ Compute the Bessel distributed amplitude within the circle, 
for m = 1 :r+l ; 

for n * l:r+l; 

x = sqrt((m-l)‘2 + (n-l)‘2); 
if x <= r; 

temp(m,n)=besseln(0,a*x) ; 
end; 
end; 

end; 

'/, Generate the entire base x base input array. 
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Y(NO:base,NO:base) = temp; 
Y(2:N0,N0:base) = f lipud(temp) ; 
Y(2:N0,2:N0) = rot90(temp ,2) ; 
Y(N0:base,2:N0) = fliplr(temp) ; 

'/, To test and view the input function: 



mesh(Y) 
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Appendix E 

Examples of Dicer 
Representations of Output Data 



This appendix contains more representations of the data from the Spyglass Dicer program. 
The representations are discussed in the text on page 4.3. 

Figures E.l and E.2 show representations of data from a Table excitation that is 23 
samples wide. Figures E.3 and E.4 show representations of data from a Circle excitation 
that is 23 samples wide. Figures E.5 and E.6 show representations of data from a Circle 
excitation that is 31 samples wide. 
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Figure E.l: Table output for d = 23 samples. Dicer representation 
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Figure E.2: Table output for d = 23 samples. Alternative Dicer representation. 
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Figure E.3: Circle output for d = 23 samples. Dicer representation. 
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Figure E.4: Circle output for d = 23 samples. Alternative Dicer representation. 
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Figure E.5: Circle output for d = 31 samples. Dicer representation. 
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No 




-0.1147 



amplitude 



2.145 



Propagation of drde(128,31) 



Figure E.6: Circle output for d = 31 samples. Alternative Dicer representation 
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